Bootstrapped outlier detection for CTD-χpod¶

What is “GE”?

Pick a bin size
In a bin we throw out top 0.6% of estimates.
How do we determine that 0.6% threshold?
1. If we assume χ is lognormally distributed, this corresponds to
2. Assume sample mean, sample std are good estimate of the population mean, population std
3. Anything greater than (sample mean) + 2.698 * (sample std) is an outlier

A bootstrap-based parametric way would be

For each resampling of the data, find which threshold value marks the top 0.7% estimate i.e. 0.7% of the data is above this threshold.
Do N resamplings so we get N estimates of the threshold.
Somehow collapse those N estimates to one value that we use to flag the data.

We can also bootstrap the GE thing by using bootstrap estimates of mean + 2.698*std

import hvplot.xarray
import matplotlib as mpl
import matplotlib.pyplot as plt
import numpy as np
import tqdm.notebook
from arch.bootstrap import IIDBootstrap
from scipy.stats import bootstrap, distributions, lognorm, norm

plt.rcParams["figure.dpi"] = 120
plt.style.use("seaborn")

Idealized dataset¶

Generate some data¶

Assume 50m bin at 0.5m resolution

N = int(50 / 0.5)
Nbad = 5

rng = np.random.default_rng()
dist = lognorm(s=2, loc=0, scale=np.exp(1e-10))
data = 1e-10 + 1e-7 * dist.rvs(size=N, random_state=rng)

# assign some outliers
idxbad = np.random.randint(0, N - 1, Nbad)
data[idxbad] = np.arange(1, Nbad + 1) * 1e-4

plt.hist(np.log10(data))

(array([ 1.,  9., 18., 25., 28., 13.,  1.,  1.,  1.,  3.]),
 array([-9.17861925, -8.59086032, -8.0031014 , -7.41534247, -6.82758355,
        -6.23982462, -5.6520657 , -5.06430677, -4.47654785, -3.88878892,
        -3.30103   ]),
 <BarContainer object of 10 artists>)

_images/bootstrap-ctd-chipod-estimates_4_1.png

Use the cdf to check the 2.698 threhsold. Looks like 0.6%

(1 - dist.cdf(dist.mean() + 2.698 * dist.std())) * 100

0.5929244418601542

Write a function that chooses outliers¶

def topk_thresh(array, thresh=0.6):
    """
    Picks value corresponding to top `thresh`% of samples.

    Effectively this is the lowest of the top `thresh`% of values.
    """

    array = array[~np.isnan(array)]
    N = len(array)
    Nthresh = int(np.ceil(thresh / 100 * N)) + 1

    if Nthresh == 1:
        raise ValueError(
            f"Looking for top N=1 point, threshold of {thresh*100}% seems low."
        )

    # sort the data and pick Nthresh-point from the high end
    return np.sort(array)[-Nthresh]


def ge_lognorm(array):
    """Assumes lognormal distribution and computes mean + 2.698*std as threshold."""
    array = array[~np.isnan(array)]
    return array.mean() + 2.698 * array.std()

Bootstrap samples of the threshold¶

v1: Use scipy bootstrap¶

This one calculates a confidence interval but that doesn’t seem so useful.

bootstrap(
    (data,),
    topk_thresh,
    n_resamples=1000,
    vectorized=False,
)

BootstrapResult(confidence_interval=ConfidenceInterval(low=2.8039623815065367e-06, high=0.0005), standard_error=0.00016452270006917774)

v2: Just generate samples of the threshold¶

This way we apply the function topk_thresh to reps number of resamplings of data. This yields reps estimates of the threshold

boot = IIDBootstrap(data)
thresh_estimates = boot.apply(topk_thresh, reps=1000).squeeze()
with np.printoptions(threshold=100):
    print(thresh_estimates)

[0.0004 0.0005 0.0005 ... 0.0005 0.0005 0.0002]

Histogram of those threshold estimates

plt.hist(np.log10(thresh_estimates))
plt.gcf().set_size_inches((6, 3))

_images/bootstrap-ctd-chipod-estimates_15_0.png

Final threshold¶

We need to pick one value from all bootstrapped estimates. Given 1000 estimates of the threshold lets pick one

median = np.median(thresh_estimates)
median

0.0004

top5 = np.sort(thresh_estimates)[-int(np.ceil(5 / 100 * len(thresh_estimates)))]
top5

0.0005

top20 = np.sort(thresh_estimates)[-int(np.ceil(20 / 100 * len(thresh_estimates)))]
top20

0.0005

ge_thresh = data.mean() + 2.796 * data.std()
ge_thresh

0.00019887193406118779

plt.hist(np.log10(data), label="log10(data)", histtype="step", lw=2)
for line, color, label in zip(
    [ge_thresh, median, top5, top20],
    ["C0", "C1", "C2", "C3"],
    ["GE", "median boot", "top5 boot", "top20 boot"],
):
    plt.axvline(np.log10(line), color=color, label=label)
plt.legend(bbox_to_anchor=(1, 1))
plt.gcf().set_size_inches((6, 3))

_images/bootstrap-ctd-chipod-estimates_21_0.png

Attempt 1¶

binned = (
    profiles.reset_coords("pres")
    .rolling(z=80, center=True)
    .construct("bin")
    .set_coords("pres")
)

OK OK these are not uniformly spaced but close enough?

(binned.pres.max("bin") - binned.pres.min("bin")).plot()

process_segment(binned.isel(z=1250).chi, debug=True)

Eddy Diffusivity

Bootstrapped outlier detection for CTD-χpod

Contents

Bootstrapped outlier detection for CTD-χpod¶

Idealized dataset¶

Generate some data¶

Write a function that chooses outliers¶

Bootstrap samples of the threshold¶

v1: Use scipy bootstrap¶

v2: Just generate samples of the threshold¶

Final threshold¶

CTD-χpod data¶

Read data¶

Process a single profile¶

Study thresholds for a single profile¶

Attempt 1¶